Использовалась библиотека stats
library(Rssa)
library(ggplot2)
library(mFilter)
data(USUnemployment)
к 14.03
Растекание — эффект, возникаемый, если длинна ряда не будет кратна периоду. Найдем длинну нашего ряда:
w <- 0.1 # Тут задаем частоту, которая соответствует периоду T = 10
x <- 1:100
y <- sin(2 * pi * w * x)
length(y)
## [1] 100
Получаем длинну ряда 100, построим периодограмму для ряда у которого длинна кратна периоду:
sp <- spec.pgram(y, spans = NULL, detrend = FALSE, log = "no", fast = FALSE, taper = 0)
Видим, что эффекта растекания на периодограмме нет, пик единственный и соответсвует заданной частоте. Теперь попробуем изменить длинну исходного ряда. Можно добавить несколько нулевых наблюдений в ряд, но лучше убрать уже имеющиеся, чтобы не искажать картину нулями:
yCut <- y[1:(length(y) - 3)] # Убираем три наблюдения
length(yCut)
## [1] 97
Теперь длинна ряда 97, что не кратно периоду T = 10, теперь посмотрим на его периодограмму:
sp <- spec.pgram(yCut, log = "no", detrend = FALSE, fast = FALSE, taper = 0)
Видим, что пик начал расползаться по соседним частотам.
Шум — это случайная составляющая временного ряда. Если средний вклад всех частот одинаковый, то шум называется белым. Свойство периодограммы белого шума: значение периодограммы имеет экспоненциальное распределение с одним и тем же средним.
n <- 1000
wnoise <- rnorm(n, 0, 1)
spec.pgram(wnoise, detrend = FALSE, log = "no", fast = FALSE, pad = FALSE,taper = 0)
Так как для шума вообще нет детерминирующих частот, мы видим из периодограммы, что каждая частота входит примерно с одинаковой силой. С увеличением количества наблюдений будет увеличиваться значение каждой частоты на сетке:
nLarge <- 10000
wnoiseLarge <- rnorm(nLarge, 0, 1)
spec.pgram(wnoiseLarge, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE, taper = 0)
Отличие красного шума от белого в том, что предыдущие значения влияют на следующие, т.е. имеют ненулевую корреляцию, в отличие от белого шума.
\(r_n = cor \cdot r_{n-1}+\sqrt{(1-cor^2)} \cdot \omega\)
w0 <- wnoise[1]
wnoise <- wnoise[2:n]
cor <- 0.8
rnoise = Reduce(function(prev_v, next_v) cor * prev_v + next_v * sqrt(1 - cor^2), wnoise, w0, accumulate = T)
spec.pgram(rnoise, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)
О данных: https://link.springer.com/chapter/10.1007/978-1-4612-5098-2_66 408 ежемесячных наблюдений безработицы в США среди мужчин Длину окна берем равную 12 в соотвествии с годовой периодичностью. Посмотрим на их график:
library(Rssa)
data(USUnemployment)
unempl.male <- USUnemployment[, "MALE"]
plot(unempl.male, type = "l")
По графику можно сказать, что временной ряд достаточно волатильный, имеет сложный тренд, по периодограмме можно сделать выводы о частотах входящих в него:
spec.pgram(unempl.male, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)
Видим, что наибольший вклад имеют низкие частоты, которые мы относим к шуму, также есть пики в точках w = 1/12 и w = 2/12, что говорит нам о наличии годовой и полугодовой сезонности.
к 14.03
movingAverage <- function(x, n=1, centered=TRUE) {
if (centered) {
before <- floor ((n-1)/2)
after <- ceiling((n-1)/2)
} else {
before <- n-1
after <- 0
}
s <- rep(0, length(x))
count <- rep(0, length(x))
new <- x
count <- count + !is.na(new)
new[is.na(new)] <- 0
s <- s + new
i <- 1
while (i <= before) {
new <- c(rep(NA, i), x[1:(length(x)-i)])
count <- count + !is.na(new)
new[is.na(new)] <- 0
s <- s + new
i <- i+1
}
i <- 1
while (i <= after) {
new <- c(x[(i+1):length(x)], rep(NA, i))
count <- count + !is.na(new)
new[is.na(new)] <- 0
s <- s + new
i <- i+1
}
s/count
}
afc <- function(filter, omega)
{
k <- seq_along(filter) - 1
h <- function(o) sum(rev(filter)*exp(-k*1i*o))
abs(sapply(omega, h))
}
n = 100
x <- 1:n
period <- 10
modelTs <- (2*x + 3) +
5*cos(2*pi*(1/period)*x) + 5*cos(2*pi*(2/period)*x) + 5*cos(2*pi*(3/period)*x) + 5*cos(2*pi*(4/period)*x) +
+ rnorm(n)
plot(modelTs, type = "l", main = "Исходный временной ряд")
sp <- spec.pgram(modelTs, detrend = FALSE, log = "no", fast = FALSE, pad = FALSE, taper = 0)
plot(modelTs, main = "Исходный ряд + Cкользящее среднее", type = "l")
legend(x = "topleft", # Position
legend = c("timeseries", "MA(5)", "MA(10)"), # Legend texts
lty = c(1, 1, 1), # Line types
col = c("black", "red", "blue"), # Line colors
lwd = 2) # Line width
maHalf <- movingAverage(modelTs, 5, TRUE)
lines(maHalf, col = "red")
maFull <- movingAverage(modelTs, 10, TRUE)
lines(maFull, col = "blue")
sp <- spec.pgram(maHalf ,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)
sp <- spec.pgram(maHalf ,detrend = FALSE, log = "yes",fast = FALSE, pad = FALSE,taper = 0)
sp <- spec.pgram(maFull ,detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)
sp <- spec.pgram(maFull ,detrend = FALSE, log = "yes",fast = FALSE, pad = FALSE,taper = 0)
freq <- seq(0, pi, 0.001)
filt <- rep(1/period, period)
omega <- freq/2/pi
str <- c("0", "1/10", "2/10", "3/10", "4/10","5/10")
plot(afc(filt, freq) ~ omega, type = "l", xlab = "Frequency", ylab = "Frequency response", xaxt = "n")
title(main = "АЧХ фильтра скользящего среднего MA(10)")
axis(1, at = seq(0, 0.5, by = 1/10), las = 2, labels = sprintf("%s", str))
period = 5
filt <- rep(1/period, period)
plot(afc(filt, freq) ~ omega, type = "l", xlab = "Frequency", ylab = "Frequency response")
title(main = "АЧХ фильтра скользящего среднего MA(5)")
График безработицы в США Периодограмма безработицы в США
Для выделения тренда выставим длинну окна равную периоду ряда (T = 12). Наложим график скользящего среднего на исходный:
plot(unempl.male, main = "Исходный ряд + Cкользящее среднее", type = "l")
legend(x = "topleft",
legend = c("timeseries", "MA(12)", "MA(24)"),
lty = c(1, 1, 1, 1),
col = c("black", "blue", "red"),
lwd = 2)
ma12 <- movingAverage(unempl.male, 12, TRUE)
ma24 <- movingAverage(unempl.male, 24, TRUE)
lines(ma12, col = "blue")
lines(ma24, col = "red")
Теперь построим периодограмму скользящего среднего:
spec.pgram(ma12, detrend = FALSE, log = "no",fast = FALSE, pad = FALSE,taper = 0)
к 14.03
Применим SSA к реальным данным, ширину окна можно взять L = K/2 и кратную периоду (условие ассимптотической разделимости). Построим матрицу взвешенных корреляций, чтобы проанализировать компоненты:
tsSSA <- ssa(unempl.male, L = length(unempl.male) %/% 2)
tsSSACor <- wcor(tsSSA, groups = 1:30)
plot(tsSSACor)
Видим, что после 4 компоненты начинается зашумленность, при этом 15-я компонента достаточно хорошо отделяется от соседних.
plot(tsSSA, type="vectors", idx = 1:20)
По парам собственных векторов можно сказать о векторах периодичностей:
plot(tsSSA, type="paired", idx = 2:20)
Видим правильные фигуры в компонентах (5, 6) и (12, 13):
tsSSASeason <- reconstruct(tsSSA, groups = list(season = c(5, 6, 12, 13)))
spec.pgram(tsSSASeason$season, log = 'no', fast = FALSE, taper = 0, detrend = FALSE)
Видим, что выделились пики соответсвующие годовой и полугодовой периодичности, как на исходной периодограмме.
tsSSASignal <- reconstruct(tsSSA, groups = list(c(1:4, 7:10), c(5, 6, 12, 13)))
plot(tsSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))
Алгоритм улучшает сильную разделимость.
tsFOSSA <- fossa(tsSSA, nested.groups = c(1:13))
tsFOSSACor <- wcor(tsFOSSA, groups = 1:30)
plot(tsFOSSACor)
plot(tsFOSSA, type = "vectors", idx=1:15)
tsFOSSASignal = reconstruct(tsFOSSA, groups <- list(c(5:13), 1:4))
plot(tsFOSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))
plot(tsSSACor)
plot(tsFOSSACor)
matplot(data.frame(tsSSASignal$F2, tsFOSSASignal$F2), type = 'l', col=c("red","green"), lty=c(1,1))
matplot(data.frame(c(unempl.male), tsSSASignal$F1, tsFOSSASignal$F1), type = 'l', col=c("black", "red","green"), lty=c(1,1))
к 21.03, 28.03
trend.poly3 = lm(unempl.male ~ poly(1:408, degree = 3))
trend.poly5 = lm(unempl.male ~ poly(1:408, degree = 5))
trend.poly7 = lm(unempl.male ~ poly(1:408, degree = 7))
plot(c(unempl.male), main = "Исходный ряд + POLY", type = "l")
legend(x = "topleft",
legend = c("timeseries", "POLY(3)", "POLY(5)", "POLY(7)"),
lty = c(1, 1, 1, 1),
col = c("black", "red", "blue", "green"),
lwd = 2)
lines(trend.poly3$fitted.values, col="red", type = "l")
lines(trend.poly5$fitted.values, col="blue", type = "l")
lines(trend.poly7$fitted.values, col="green", type = "l")
LOESS метод основан на сглаживании с помощью построения локальных взвешенных линейных регрессий
index <- 1:408
fit.loess <- loess(unempl.male ~ index, span = 0.2, degree = 1)
trend.loess <- predict(fit.loess)
plot(c(unempl.male), main = "Исходный ряд + LOESS", type = "l")
legend(x = "topleft",
legend = c("timeseries", "LOESS"),
lty = c(1, 1, 1, 1),
col = c("black", "blue"),
lwd = 2)
lines(trend.loess, col = "blue", type = "l")
chest = hpfilter(unempl.male, freq=60000, type="lambda")
trend.hp = chest$trend
plot(unempl.male, type="l")
lines(trend.hp, col="red")
legend("topleft", legend = c("Real TS","HP"), col = c("black","red"), lty=1)
plot(unempl.male, type="l")
lines(trend.hp, col="red")
lines(ma24, col="blue")
lines(ts(trend.loess, start=c(1948, 1), frequency=12), col="orange")
lines(ts(trend.poly7$fitted.values, start=c(1948, 1), frequency=12), col="green")
legend("topleft", legend = c("Real TS","HP", "MA(24)", "LOESS", "POLY(5)"), col = c("black","red", "blue", "orange", "green"), lty=1)
04.04
STL имеет следующие параметры:
(inner) – число итераций внутреннего цикла (р).
(outer) – число итераций внешнего цикла. У нас в данных нет аутлаеров, поэтому outer = 0
(l.window) – сглаживающий параметр для low-pass фильтра. Возьмем равным 13. Рекомендуют брать ближайший нечетный к периоду.
(t.window) – сглаживающий параметр для тренда. Рекомендуют брать ближайшим нечетным к (1.5*period) / (1-(1.5/s.window)). Равен 23.
ns - (s.window) – сглаживающий параметр сезонности. Берем 13, так как годовая периодичность и параметр должен быть нечетным
tsSTL <- stl(unempl.male, s.window = 13, l.window = 13, outer = 0, inner = 1, t.window = 23)
plot(tsSTL)
tsCSD <- ssa(unempl.male, L = length(unempl.male) %/% 2, force.decompose = FALSE, svd.method = "nutrlan")
fit.dec <- decompose(tsCSD, neig = 40)
tsCSDCor <- wcor(fit.dec, groups= 1:40)
plot(tsCSDCor)
plot(fit.dec, type = "vectors", idx = 1:20)
tsCSDSignal <- reconstruct(fit.dec, groups = list(trend = c(1:4, 7:10, 15), season = c(5, 6, 12, 13)))
plot(tsCSDSignal)
plot(unempl.male, type="l")
lines(tsSSASignal$F1, col="red")
lines(tsFOSSASignal$F1, col="blue")
lines(tsCSDSignal$trend, col="orange")
lines(ts(data.frame(tsSTL$time.series)$trend, start=c(1948, 1), frequency = 12), col="green")
legend("topleft", legend = c("Real TS","SSA", "FOSSA", "CSD", "STL"), col = c("black","red", "blue", "orange", "green"), lty=1)
Projection SSA применяется для улучшения разделимости линейного тренда и сезонности. Используем смоделированные данные с линейным трендом (Применение фильтра скользящего среднего к модельным данным):
x <- 1:100
modelTs <- (2*x + 3) +
+ 5*cos(2*pi*(1/25)*x) + cos(2*pi*(1/5)*x) + rnorm(50)
plot(modelTs, type="l")
# Применим SSA с двойным центрированием для лучшего отделения линейного тренда
modelTsSSA <- ssa(modelTs, L = 25, column.projector='centering', row.projector='centering')
plot(modelTsSSA, type = "vectors")
plot(modelTsSSA, type = "paired")
modelTsSSASignal <- reconstruct(modelTsSSA, groups = list(1:2))
plot(modelTs, type = "l")
lines(modelTsSSASignal$F1, col = "red")
Теплицев SSA применяется как улучшение basic-SSA для стационарных рядов.
x <- 1:100
modelTs <- 5*cos(2*pi*(1/25)*x) + 2*cos(2*pi*(1/5)*x) + rnorm(50)
plot(modelTs, type="l")
modelTsSSA <- ssa(modelTs, kind = "toeplitz-ssa")
plot(modelTsSSA, type = "vectors")
plot(modelTsSSA, type = "paired")
modelTsSSASignal <- reconstruct(modelTsSSA, groups = list(1:4))
plot(modelTsSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))
plot(modelTs, type = "l")
lines(modelTsSSASignal$F1, col = "red")
02.05
forcastR <- rforecast(tsSSA, groups = list(trend = c(1:4, 7:10), season = c(5, 6, 12, 13)), len = 12, only.new = F)
plot(unempl.male, type="l")
lines(data.frame(forcastR)$trend, col="red")
legend("topleft", legend = c("Real TS","RForcast trend"), col = c("black","red"), lty=1)
Применяем SSA к ряду с отсеченным периодом
tsCut = ts(unempl.male[1:(length(unempl.male)-12)], start=c(1948, 1), frequency=12)
tsCutSSA <- ssa(tsCut, L = length(tsCut) %/% 2)
plot(wcor(tsCutSSA, groups = 1:30))
plot(wcor(tsCutSSA, groups = 1:15))
plot(tsCutSSA, type="vectors", idx = 1:20)
tsCutSSASignal = reconstruct(tsCutSSA, groups <- list(c(1:3, 6:10, 14:15), c(4,5,12,13)))
plot(tsCutSSASignal, add.residuals = TRUE, add.original = TRUE, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 2))
forcastRCut <- rforecast(tsCutSSA, groups = list(trend = c(1:3, 6:10, 14:15), season = c(4,5,12,13)), len = 12, only.new = F)
plot(c(tsCut), type="l")
lines(c(data.frame(forcastRCut)$trend), col="red")
lines(c(data.frame(forcastR)$trend), col="blue")
legend("topleft", legend = c("Cut TS","RForcast cut trend", "RForcast full trend"), col = c("black","red", "blue"), lty=1)
forcastV <- vforecast(tsSSA, groups = list(trend = c(1:4, 7:10), season = c(5, 6, 12, 13)), len = 12, only.new = F)
plot(cbind(unempl.male, forcastV$trend), plot.type = "single",
col = c("black", "red"), ylab = NULL)
legend("topleft", legend = c("Real TS","VForcast full trend"), col = c("black","red"), lty=1)
forcastVCut <- vforecast(tsCutSSA, groups = list(trend = c(1:3, 6:10, 14:15), season = c(4,5,12,13)), len = 12, only.new = F)
plot(c(tsCut), type="l")
lines(c(data.frame(forcastVCut)$trend), col="red")
lines(c(data.frame(forcastV)$trend), col="blue")
legend("topleft", legend = c("Real TS","VForcast cut trend", "VForcast full trend"), col = c("black","red", "blue"), lty=1)
25.04/02.05
library(Matrix)
Находим ранг матрицы для экспоненциального ряда:
N = 10
LRF = 3.1^(1:N)
rankMatrix(hankel(LRF))[1]
## [1] 1
Ожидаемо, он равен одному. Далее находим сигнальные корни ряда:
LRF.ssa = ssa(LRF, L=2, method = "svd")
l = lrr(LRF.ssa, groups = list(1))
plot(l)
roots_my = roots(l)
K = length(roots_my)
S_n = LRF[1:K]
lin_sys = data.frame(S_n)
vars = matrix(nrow = K, ncol = K)
for(i in 1:K){
vars[i]=roots_my[i]^(1:K)
}
lm0 = lm(formula = S_n ~ 0 + ., data = data.frame(vars))
print(lm0)
##
## Call:
## lm(formula = S_n ~ 0 + ., data = data.frame(vars))
##
## Coefficients:
## vars
## 1
Находим значения по линейной рекурентной формуле:
pred_S_n = lm0$coefficients[[1]]*(roots_my[1]^(1:N))
print(LRF)
## [1] 3.1000 9.6100 29.7910 92.3521 286.2915 887.5037
## [7] 2751.2614 8528.9104 26439.6222 81962.8287
S_n = LRF[1:N]
lin_sys = data.frame(S_n)
vars = matrix(nrow = N, ncol = K)
for(i in 1:K){
vars[i]=roots_my[i]^(1:N)
}
lm0 = lm(formula = S_n ~ 0 + ., data = data.frame(vars))
print(lm0)
##
## Call:
## lm(formula = S_n ~ 0 + ., data = data.frame(vars))
##
## Coefficients:
## vars
## 1
Сравниваем исходный ряд с полученным
pred_S_n = lm0$coefficients[[1]]*(roots_my[1]^(1:N))
print(LRF)
## [1] 3.1000 9.6100 29.7910 92.3521 286.2915 887.5037
## [7] 2751.2614 8528.9104 26439.6222 81962.8287
print(pred_S_n)
## [1] 3.1000 9.6100 29.7910 92.3521 286.2915 887.5037
## [7] 2751.2614 8528.9104 26439.6222 81962.8287
N = 10
LRF = 2*(1:N)+5
Находим ранг матрицы:
rankMatrix(hankel(LRF))[1]
## [1] 2
Линейная функция задается двумя компонентами.
LRF.ssa = ssa(LRF, L=3, method = "svd")
l <- lrr(LRF.ssa, groups = list(1:2))
roots_my <- roots(l)
print(l)
## [1] -1 2
## attr(,"class")
## [1] "lrr"
print(roots_my)
## [1] 1.0000002 0.9999998
Получаем два комплексно-сопряженных корня
print(2 * pi / Arg(roots_my))
## [1] Inf Inf
print(Mod(roots_my))
## [1] 1.0000002 0.9999998
roots_my2 <- rep(Re(mean(roots_my)), length(roots_my))
print(roots_my2)
## [1] 1 1
Находим коэффиценты исходного ряда
m = length(roots_my2)
s_n = LRF[1:N]
vars <- matrix(nrow = N, ncol = m)
for (i in 1:m) {
vars[ ,i] <- (1:N) ^ (i-1) * roots_my2 ^ (1:N)
}
lm0 <- lm(s_n ~ 0 + ., data = data.frame(vars))
print(lm0)
##
## Call:
## lm(formula = s_n ~ 0 + ., data = data.frame(vars))
##
## Coefficients:
## X1 X2
## 5 2
Сравниваем исходный ряд с полученным
pred_S_n = lm0$coefficients[[1]][1] * (vars[,1]) + lm0$coefficients[[2]][1] * (vars[,2])
print(LRF)
## [1] 7 9 11 13 15 17 19 21 23 25
print(pred_S_n)
## [1] 7 9 11 13 15 17 19 21 23 25
Проанализируем собственные числа ряда
tsSSA$sigma
## [1] 367075.857 74505.552 63304.591 35713.403 34728.609 34583.735
## [7] 27721.976 22861.296 22230.219 18781.365 18634.256 18365.270
## [13] 18171.464 14245.491 13430.646 13113.963 10353.466 10239.538
## [19] 9456.080 9280.711 9009.749 8960.235 8234.105 8146.040
## [25] 7272.790 6845.598 6199.300 6024.433 5492.880 5475.516
## [31] 5071.906 4194.408 4169.282 3688.478 3611.481 3376.238
## [37] 3365.972 3203.616 3192.859 3179.853 3072.833 3014.386
## [43] 2996.338 2987.268 2984.135 2906.750 2672.857 2504.003
## [49] 2419.366 2386.851
plot(tsSSA$sigma, type="l")
Видим 13 “ненулевых” собственных числа.
rk = 13
par = parestimate(tsSSA, groups = list(1:rk), method = "esprit")
modulusReal = par$moduli
periodsReal = par$periods
o = order(abs(periodsReal), decreasing = TRUE)
r0 = reconstruct(tsSSA, groups = list(signal = 1:rk))$signal
len <- rk
vars <- matrix(nrow = len, ncol = rk)
for (i in 1:rk) {
if (periodsReal[i] == Inf)
vars[, i] <- modulusReal[i]^(1:len)
else if (periodsReal[i] == 2)
vars[, i] <- (-modulusReal[i])^(1:len)
else if (periodsReal[i] > 0)
vars[, i] <-
modulusReal[i]^(1:len) * sin(2 * pi * (1:len) / periodsReal[i])
else
vars[, i] <-
modulusReal[i]^(1:len) * cos(2 * pi * (1:len) / periodsReal[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
print(round(coefs0), digits = rk)
## X1 X2 X3 X4 X5 X6 X7 X8
## -4965469 -5555248 7213618 54636156 -48963499 1651 739 205
## X9 X10 X11 X12 X13
## -65 1225627 -117262 NA NA
idx <- seq(1, rk)
coefs.c.phase <- numeric(length(idx))
phases.c <- numeric(length(idx))
periods.c.phase <- numeric(length(idx))
moduli.c.phase <- numeric(length(idx))
for (i in seq_along(idx)) {
periods.c.phase[i] <- periodsReal[idx[i]]
moduli.c.phase[i] <- modulusReal[idx[i]]
coefs.c.phase[i] <- coefs0[idx[i]]
phases.c[i] <- atan2(coefs0[idx[i] + 1], coefs0[idx[i]])
if(i == rk){
phases.c[i] <- atan2(coefs0[idx[i]], coefs0[idx[i]])
}
}
s = data.frame(periods = periods.c.phase, phases = phases.c,
coefficients = coefs.c.phase,
moduli =moduli.c.phase)
print(s)
## periods phases coefficients moduli
## 1 62.059089 -2.30019430 -4.965469e+06 1.0102164
## 2 -62.059089 2.22704139 -5.555248e+06 1.0102164
## 3 207.578092 1.43952547 7.213618e+06 1.0031497
## 4 -207.578092 -0.73069723 5.463616e+07 1.0031497
## 5 Inf 3.14155894 -4.896350e+07 1.0007054
## 6 11.995011 0.42101827 1.650952e+03 0.9986643
## 7 -11.995011 0.27033942 7.392869e+02 0.9986643
## 8 6.007629 -0.30780195 2.048739e+02 0.9982437
## 9 -6.007629 1.57084947 -6.513055e+01 0.9982437
## 10 38.009961 -0.09538464 1.225627e+06 0.9969873
## 11 -38.009961 NA -1.172619e+05 0.9969873
## 12 50.401313 NA NA 0.9907013
## 13 -50.401313 NA NA 0.9907013
ts1 <- read.csv(file = "ts1.txt", header = TRUE, as.is = FALSE)
ts5 <- read.csv(file = "ts5.txt", header = TRUE, as.is = FALSE)
plot(ts1$x, type = "l")
plot(ts5$x, type = "l")
acf(ts1, main = "Функция для данных ts1")
acf(ts5, main = "Функция для данных ts5")
По количеству столбцов, значимо отличающихся от нуля, можно определить количество параметров в авторегрессии.
pacf(ts1, main = "Функция для данных ts1")
pacf(ts5, main = "Функция для данных ts5")
fit1 <- auto.arima(ts1)
fit5 <- auto.arima(ts5)
fit1
## Series: ts1
## ARIMA(2,0,2) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 mean
## -0.2652 0.6067 1.1344 0.1612 99.9179
## s.e. 0.0481 0.0335 0.0542 0.0454 0.1132
##
## sigma^2 = 1.064: log likelihood = -1448.07
## AIC=2908.15 AICc=2908.23 BIC=2937.59
fit5
## Series: ts5
## ARIMA(3,0,0) with zero mean
##
## Coefficients:
## ar1 ar2 ar3
## 0.3103 -0.1605 0.0984
## s.e. 0.0222 0.0230 0.0223
##
## sigma^2 = 0.3961: log likelihood = -1910.39
## AIC=3828.78 AICc=3828.8 BIC=3851.18
forecastedValues1 <- forecast(fit1, 20)
forecastedValues5 <- forecast(fit5, 40)
plot(forecastedValues1, main = "Graph with forecasting FIT1",
col.main = "darkgreen")
plot(forecastedValues5, main = "Graph with forecasting FIT5",
col.main = "darkgreen")
fitArima <- auto.arima(unempl.male)
checkresiduals(fitArima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2)(2,1,1)[12] with drift
## Q* = 20.36, df = 17, p-value = 0.2562
##
## Model df: 7. Total lags used: 24
summary(fitArima)
## Series: unempl.male
## ARIMA(2,0,2)(2,1,1)[12] with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sar2 sma1 drift
## 1.6149 -0.6340 -0.4319 0.1149 -0.0108 0.1072 -0.8538 5.8252
## s.e. 0.1395 0.1366 0.1483 0.0671 0.1022 0.0932 0.0909 3.6365
##
## sigma^2 = 16160: log likelihood = -2484.99
## AIC=4987.98 AICc=4988.44 BIC=5023.81
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.6993681 123.9684 91.59157 -0.34855 5.249927 0.1989636
## ACF1
## Training set 0.0005744349
fitArimaForecast <- forecast(fitArima, h = 24)
autoplot(fitArimaForecast)
tsSSABforecast <- bforecast(tsSSA, groups = list(c(1:4, 7:10)), len = 24, R = 50, only.new = F)
tsSSAETS <- hw(ts(unempl.male[1:(length(unempl.male)-24)], start=c(1948, 1), frequency=12), seasonal="multiplicative", h = 24)
autoplot(tsSSABforecast)
plot(tsSSAETS)
sd(unempl.male)
## [1] 794.9778
tsTest = tail(unempl.male, 24)
tsSSABforecast <- bforecast(tsSSA, groups = list(c(1:4, 7:10)), len = 24, R = 50, only.new = T)
mseSSA = mean((data.frame(tsSSABforecast)$Value - c(tsTest))**2)
mseARIMA = mean((c(fitArimaForecast$mean) - c(tsTest))**2)
mseETS = mean((c(tsSSAETS$mean) - c(tsTest))**2)
print(sqrt(mseSSA))
## [1] 1399.446
print(sqrt(mseARIMA))
## [1] 1024.97
print(sqrt(mseETS))
## [1] 889.6957
library(Rssa)
library("lattice")
firstHalf = c(10*sin(2*pi*(1:72)/6))
secondHalf = c(30*sin(2*pi*(73:144)/6))
tsdestruct = c(firstHalf,secondHalf)
plot(tsdestruct, type = "l", col = "red")
s <- ssa(tsdestruct, L = 12)
w <- wcor(s, groups = 1:10)
plot(w)
plot(s, type = "vectors", idx=1:10)
r <- reconstruct(s, groups = list(Trend = c(1)))
tsdestruct_res <- residuals(r)
N <- length(tsdestruct_res)
rank <- 2
periods <- function(M, L) {
ts(sapply(1:(N - M),
function (i) {
s <- ssa(tsdestruct_res[i:(i + M - 1)], L = L)
par <- parestimate(s, groups = list(c(1:rank)),
method = "esprit")
abs(par$periods[1])
}),
start = time(tsdestruct)[M + 1], delta = 1)
}
per12 <- periods(12, 6)
per24 <- periods(24, 12)
lattice::xyplot(plot.method = "xyplot",per12, type = "l")
M <- 12; L <- M / 2
hm <- hmatr(tsdestruct_res, B = M, T = M, L = L, neig = rank)
plot(hm)